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We reconstruct the trajectories of ultra-high energy cosmic rays (UHECR) — observed by the 
AGASA experiment — in the Galactic magnetic field assuming that all particles have the same 
charge. We then study correlations between the reconstructed events and BL Lacs. The correlations 
have significance below 10 -3 in the case of particles with charge +1. In the case of charge —1 the 
correlations are absent. We interpret this as evidence that protons are present in the flux of UHECR. 
Observed correlation provides an independent evidence that BL Lacs emit UHECR. 



I. INTRODUCTION 

Ultra-high energy cosmic rays (UHECR) are a subject 
of an active debate for over more than 20 years. The data 
accumulated during this time f|, ^ provide a compelling 
evidence that the GZK cutoff 0] predicted in the UHECR 
spectrum at energies E ~ 10 eV may be absent. Many 
models explaining this puzzle have been proposed (for 
reviews see Refs. Q). The actual focus of the debate is 
the question of whether an unconventional astrophysical 
model can be constructed which explains the observed 
super-GZK events, or a new physics is required. This 
key question still remains open. 

With the accumulation of the high-energy events im- 
portant signatures emerge which discriminate efficiently 
between different models. It has been known for some 
time that UHECR events form clusters ||, [| Recent 
analysis shows that clustering is statistically significant 
0. The small angular size of clusters of order ~ 2.5°, 
consistent with the experimental angular resolution, sug- 
gests that they are due to point sources. The models 
which do not reproduce this feature are therefore disfa- 
vored. 

Observation of clusters implies that some sources of 
UHECR may be identified. Recently, significant correla- 
tions of arrival directions of UHECR with positions of 
BL Lacertae (BL Lacs) were found ||. BL Lacs are 
blazars (ANGs with relativistic jet directed along the line 
of sight) characterized in particular by the (near) absence 
of the emission lines. The correlations between UHECR 
and BL Lacs are most significant at angles of order ~ 2.5° 
and are present at relatively low energies E > 2.4 x 10 19 . 
Such tightness suggests strongly the existence of neutral 
primary particles. Association with distant BL Lacs com- 
bined with neutrality of primaries rules out most of the 
models of UHECR, leaving the models based on neutrino 
H (the Z-burst models), models with hypothetical "im- 
mune messengers" [[To[, or models involving violation of 
the Lorentz invariancc [[fl] . It should be noted that ultra- 
high energy photons also cannot be excluded at present 



Regardless of whether there exist neutral primary par- 
ticles, correlations with BL Lacs imply that the acceler- 
ation mechanism of UHECR production actually works. 
Therefore, protons are involved and should be present in 
the UHECR flux at energies around or below the GZK 
cutoff. If identified, such protons provide independent 
evidence that BL Lacs are indeed sources of UHECR. 
The purpose of this paper is to address this issue. 

Knowing actual sources of UHECR is an extremely 
powerful tool even when statistics is limited. As we will 
see shortly, this tool can be used to study charge com- 
position of UHECR. The idea is to use the bending of 
charged primary particles in the Galactic (GMF) |l5| and 
extragalactic (EGMF) |ll| magnetic fields. While deflec- 
tions in EGMF are random and therefore unpredictable, 
deflections in GMF are regular. If charges, energies of 
particles and GMF are known, the original directions 
(before entering GMF) can be restored. If the effect of 
EGMF is small, the actual positions of sources can be 
reconstructed. 

In practice, one has to solve the inverse problem: to re- 
construct charges assuming set of potential sources. We 
show that this is possible at least in a statistical sense: as- 
suming that a substantial fraction of UHECR are protons 
significantly improves correlations with BL Lacs. This 
implies simultaneously that GMF model is roughly cor- 
rect, and that the effect of EGMF is small. 

The paper is organized as follows. In Sect. II we briefly 
review the present knowledge about the Galactic mag- 
netic field and fix the GMF parameters for further cal- 
culations. In Sect. |H we analyse correlations between 
BL Lacs and UHECR when the latter are assigned non- 
zero charges. Sect. |IV| contains discussion and concluding 
remarks. 



II. GALACTIC MAGNETIC FIELD 

The Galactic magnetic field can be divided in two 
parts: the disc and the halo. Each one has regular and 
turbulent components. While the strength of turbulent 
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component is larger, it is the regular field which gives 
dominant contribution into CR propagation. Most of 
the information on the regular component of the disc 
is obtained from the Faraday rotation measurements of 
pulsars and extragalactic radio sources. The latter are 
used also for the reconstruction of the halo field. Mag- 
netic field in the disc resembles the spiral structure of our 
Galaxy. It may either reverse direction between different 
spiral arms (bisymmetric, or BSS model), or there may 
be no reversals (axisymmetric, or ASS model). Several 
field reversals were detected |l§ [9| [H) which 
are consistent with BSS model (note however that dis- 
crimination between the two models is complicated by 
small-scale irregularities in the magnetic field). In our 
calculations we adopt the BSS model. Simple analyti- 
cal representation of the spiral field structure (see Refs. 
||l8[ p2| ) contains the following parameters: 

• Distance from the Sun to the Galactic center, R = 
8.5 kpc. 

• Local (at the Sun position in the Galaxy) field 
strength, Bo- 
rn Pitch angle p which determines the direction of the 

local magnetic field (the field points in the direction 
I = 90 + p in the Galactic coordinates). 

• Distance d to the first field reversal. Negative d 
means that nearest reversal occurs in the direction 
to the Galactic center, positive corresponds to the 
opposite direction. 

In terms of these parameters the field in the disk is 
written as 

Bg = B cos(p), B r = B sin(p). 
The magnitude B — B(r, 9) has the spiral structure, 

B(r, 9) = B{r) cos (e - /3hi (£) + , (1) 

where (3 = l/tan(p), a constant phase <fi is given by 
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In the last expression the standard assumption that mag- 
nitude of the field decreases as r _1 in radial direction is 
made. It is also assumed that B(r) — const at r < 4 
kpc. Note that precise dependence of the disc field far 
away from the Sun position is not important for our study 
as only a small fraction of the observed UHECR passes 
through this region. Note also that the constant phase 6 
can be absorbed in another parameter r n , so that Eq. (PP) 
becomes B(r,9) — B(r) cos(6> — (3 ln(r/ro)), which is the 



parametrization used in Refs. [|18|, We find it more 
convenient to work with parameters directly related to 
the local field. 

To proceed, we need to fix the parameters Bo, p and 
d. All studies of GMF based on pulsar and extra gala ctic 
rotation measures converge on B = 1.4 ^G, see |17| , [Ts| , 
[H], and recent reviews |24|. We adopt this value. 

Pitch angle was found to be close to zero or positive 
in early publications, p = +5° Q, p = -2° Q, but 
decreased in more recent studies: p = —8° was obtained 
in Refs. g ^ while in Ref. |H) p = -15° was 
found as an average pitch angle in nearby spiral arms. 
Following reviews |23|, Q we take p = —8°. 

Field reversal was found to be at d = —0.2 — 0.3 kpc 
in Ref. Jl8| (however, the best fit value of ro in the same 
paper corresponds to d = —0.48), at d = —0.4 kpc in 
Ref. @, at d = -0.6 kpc in Refs. Q H, M- Fin ally, 
in the review |2^] d = —0.6 kpc was cited, while the 
review §| follows Ref. (l8). We take d = -0.5 kpc. 

The simplest approximation for the halo field is ob- 
tained by taking the disk field and extending it outside 
of the disk with exponentially decreasing amplitude, 



B(r, 9, z) = cxp 



B(r,9) . 



(4) 



This introduces one more parameter, the height h. Small 
values of B z = 0.2/iG found in |18 for the halo field 
support this approximation. The disc has a height of 
h = 1.5 kpc according to reviews (2^, Q and we adopt 
this value. 

In addition, the halo fields above and below the disk 
(more precisely, its parts parallel to the disk) may be 
parallel or anti-parallel. In the latter case the halo field 
may be approximated by 



B(r, 9, z) = sign(z) exp ( — — 



B(r,( 



(5) 



The first case, Eq. (|f|), corresponds to the quadrupole- 
type model which we denote BSSq, while the second case, 
Eq. (||), corresponds to the dipole-type model denoted 
as BSSd in what follows. There arc indications in favour 
of the BSS D global structure ||, |§, J|| of the halo field. 

In calculations of UHECR propagation in GMF, the 
parameters of the latter are often chosen following early 
work p3[ |. We do not use the conventions of Ref. |l4|] 
because they are ambiguous: the parameter ro used there 
does not correspond to the cited local field strength, the 
value of pitch angle does not correspond to cited direction 
of the local field and references for the assumed value of 
the halo height are not given. 

To summarize, we adopt Eqs. (|]-||) for the Galactic 
magnetic field with the following set of parameters 



B = 1.4 ^tG 
d = —0.5 kpc 



P= -» 
h = 1.5 kpc 



(6) 



and assume that this field extends to R max — 20 kpc in all 
directions. It should be stressed that these parameters 



3 



are chosen on the basis of rotation measurements and 
their choice has nothing to do with the propagation of 
UHECR. 



III. CORRELATION ANALYSIS 

We start by specifying the sets of BL Lacs and 
UHECR. Following Ref. ||, we take BL Lacs form the 
QSO catalog [^9| which contains 306 confirmed BL Lacs. 
The choice of UHECR set is motivated as follows. Apri- 
ori, best results should be achieved with the largest set 
having best angular and energy resolution. In addition, 
this should be (relatively) high-energy set, because the 
uncertainties in the deflection angle (which result from 
uncertainties in GMF and energies of particles) should 
not dominate over the angular resolution. This require- 
ment is satisfied at E > 4 x 10 19 eV. Therefore, we chose 
all published AGASA events with energy E > 4 x 10 19 eV 
This set contains 57 events. 

For the quantitative measure of correlations we use 
the probability p(S) introduced in Refs. j?], |^], which is a 
probability to have certain excess of events within the an- 
gle 8 from any of the sources (BL Lacs) . This probability 
is calculated by counting how often the excess observed 
in the real data occurs in the Monte-Carlo (MC) simula- 
tions. The MC configurations are generated as described 
in Refs. jj], 0|. The difference is that now charges are as- 
signed and arrival directions are corrected for GMF both 
in the real data and in each MC set. The energies of MC 
events are taken from the real data. Exactly the same 
treatment of the real data and MC sets is important to 
prevent appearance of artificial correlations. 

The charge assignment can be done by any algorithm 
as long as the same algorithm is used in MC simulations. 
In the simplest case one assigns equal charges to all par- 
ticles. Although correlations which are due to neutral 
particles are destroyed in this case, this effect may be 
compensated by charged particles which move closer to 
their sources. One may expect this situation when frac- 
tions of neutral and charged particles are comparable. 

As explained in Ref. ^ , there are two possible strate- 
gies to estimate the significance of correlations on the 
basis of p(8). One may impose cuts on BL Lac set and 
adjust them in such a way that correlations are maxi- 
mum. Likewise, one may look for the "best" values of 
the parameters of GMF. In this case the penalty fac- 
tor should be calculated and included in the probability. 
Since reconstruction of particle trajectories in GMF is 
rather time-consuming, this approach is not very prac- 
tical in the case at hand. Thus, we do not adjust the 
parameters of GMF, Eq. (^jj), and impose no cuts on BL 
Lac set except for a single cut in the apparent magnitude. 
Instead of penalty calculation we present explicit depen- 
dence of the probability p(8) on this cut. We draw our 
main conclusions from the behaviour of the curve, rather 
than from the value of the probability at the minimum. 

Correlations found in Ref. I| were strongest at 8 ~ 




FIG. 1: The dependence of the probability p(2.5°) on the 
cut in magnitude in BL Lac catalog. Quadrupole-type GMF 
model, Eq. M), is assumed. 



2.5°. We therefore fix this value of 8 in our calculations. 
Since 8 is not adjusted to minimize the probability, there 
is no penalty factor associated with that. 

Consider first the case of the symmetric (quadrupole) 
model BSSq. Fig. |l| shows the dependence of p(2.5°) on 
the cut in magnitude (the cut mag < 20 corresponds to 
inclusion of practically all BL Lacs). Solid and dotted 
lines correspond to the charge +1 or —1 assigned to all 
particles, respectively. We see that in both cases p(2.5°) 
has minima which are comparable in depth (although 
both not very deep). One is tempted to conclude that 
both charges may be present. This kind of situation is 
expected in the Z-burst models. One may notice, how- 
ever, that these minima are not equally significant. The 
minimum at Q = +1 is wider and corresponds to much 
higher statistics (correlations are present at the level of 
~ 1% even when all BL Lacs are included). Moreover, 
event-by-event analysis shows that 12 out of 14 events 
contributing to the minimum of probability at Q = +1 
are situated in the Northern hemisphere. This strange 
feature suggests that the field in the Southern hemisphere 
is wrong, while in the Northern one it is roughly correct. 
The obvious thing to try is the BSSd model in which the 
Southern field has different sign. 

Consider the case of the asymmetric (dipole) model 
BSSd- In our problem changing the direction of the GMF 
in the Southern hemisphere to the opposite is equivalent 
to flipping sign of charges of the events in the Southern 
hemisphere. According to previous results, this should 
increase the correlations. Indeed, the situation changes. 
Fig. H shows the dependence of the probability p(2.5°) 
on the cut in magnitude in the case of the BSSd model. 
Three cases Q = —1,0, +1 are shown. The correlations in 
the case Q — +1 have improved by almost two orders of 
magnitude as compared to the BSSq model and are now 
at the level below 10 -3 in a wide range of magnitudes. 
Even with no cuts on BL Lacs the significance of corre- 
lations at Q = +1 is below 10 -3 . On the contrary, in the 
case Q = — 1 the correlations are now absent. In the case 
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FIG. 2: The dependence of the probability p(2.5°) on the cut 
in magnitude in BL Lac catalog. Dipole-type GMF model, 
Eq. (|5|), is assumed. 




FIG. 3: The dependence of the probability p(8) on the angle 
5 (in degrees) with the cut mag < 18 in the BL Lac catalog 
and charge Q = +1. Dipole-type GMF model, Eq. (^), is 
assumed. 



Q = correlations with BL Lacs satisfying mag < 18 
are observed at the level of 2%. For completeness, Fig. || 
shows the the dependence of p(S) on the angle S at the 
cut mag < 18. 



IV. DISCUSSION AND CONCLUSIONS 

The pairs BL Lac - cosmic ray separated by less than 
2.5° are listed in Table I (imposing the cut mag < 18). 
BL Lacs appearing in this table are probable sources of 
UHECR. For each BL Lac the name, Galactic coordi- 
nates 1 and b, and redshift z (when known) are given. 
Corresponding cosmic ray energy E (in units of 10 19 eV) 
and probable charge Q are listed in the last two columns. 

We have assigned Q — to those particles which con- 
tribute to correlations at 6 — 2.5° in the case when all 
particles are assumed to be neutral. Likewise, we have 
assigned charge Q — +1 to those events which fall within 
2.5° from any of the BL Lacs when they are assumed 





Name 


1 


b 


z 


E 


Q 


1 


2EG J0432+2910 


170.52 


-12.6 


- 


5.47 


or 


+ 1 


2 










4.89 


or 


+1 


3 


RX J1838.7+4802 


76.95 


21.83 


- 


10.6 


or 


+1 


4 










4.35 


+1 




5 


RGB J0109+182 


128.82 


-44.4 


- 


21.3 


+1 




6 










5.07 


or 


+1 


7 


RX J1058.6+5628 


149.59 


54.42 


0.144 


7.76 







8 










5.35 







9 


RGB J1415+485 


91.2 


63.11 


- 


6.22 


+1 




10 


RX J0035.2+1515 


117.15 


-47.44 


- 


5.53 


or 


+1 


11 


RX J1704.8+7138 


103.09 


33.96 


- 


4.78 


+ 1 




12 


OT 465 


74.22 


31.4 


- 


4.88 


or 


+1 


13 


RX J1702.6+3115 


53.4 


35.76 


- 


4.47 


or 


+1 


14 


RX J1359.8+5911 


107.36 


55.83 


- 


4.46 







15 


RGB J0159+107 


148.75 


-48.64 


- 


4.2 


+1 




16 


1ES 1853+671 


97.74 


24.63 


0.212 


4.39 


+1 




17 


RX J1100.3+4019 


175.87 


63.56 


- 


7.21 







18 


EXO 1118.0+4228 


167.85 


66.16 


0.124 




or 


+ 1 


19 


RGB J1426+340 


57.6 


68.53 




4.97 


+1 




20 


TEX 1428+370 


63.95 


66.92 


0.564 




+ 1 




21 


B2 0804+35 


186.48 


30.35 


0.082 


4.09 


+ 1 




22 


TXS 0806+315 


190.42 


29.36 


0.22 




+ 1 





TABLE I: The list of pairs BL Lac - cosmic ray which con- 
tribute to correlations of Fig. ^| at mag < 18. 



to have charge +1 and deflection in GMF is taken into 
account. Some events satisfy both requirements; corre- 
sponding entry in Table I reads "0 or +1". 

Lines 1-8 of Table I contain four BL Lacs which corre- 
late with two cosmic rays each. These are most probable 
sources. In the lines 9-16 eight BL Lacs are listed which 
correlate with singlets. Finally, lines 17-22 contain three 
cosmic rays for which sources are ambiguous (each ray 
has two neighbouring BL Lacs within S < 2.5°). 

Examining Table I one notices two striking regularities 
(we do not attempt to assign statistical significance to 
those in the present paper). First, most objects in this 
table are X-ray selected radio loud BL Lacs. Second, the 
fraction of BL Lacs with unknown redshifts in Table I 
is much larger than in average over the whole BL Lac 
catalog. In doublet part of the Table this fraction is 6:2, 
in singlet section this fraction is 7:1. 

Absence of emission lines (more precisely, their weak- 
ness and narrow width) is a defining feature of BL Lac 
family within the general blazar class. Therefore, it is not 
surprising that redshifts of roughly half of confirmed BL 
Lacs are not known. Increased fraction of such BL Lacs 
in Table I may mean that the absence of emission lines 
is important for a blazar to became emitter of UHECR. 

As it is generally assumed, BL Lacs with unknown red- 
shift may be far away: the observational bias would then 
explain that the redshifts of such objects are unknown 
more often. Then however new physics or extreme astro- 
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physics would be necessary to explain observed correla- 
tions. 

Within the conventional framework, when propagation 
is limited by interactions with comic backgrounds, one 
should expect sources to be relatively nearby. Closest 
BL Lac with known redshift is at z — 0.029. This is well 
outside of the GZK sphere. It should be noted, however, 
that most of the cosmic rays in Table I have energies 
below the GZK cutoff. As it was found in Ref. fj^ ], the 
flux of protons at energy £ = 8x 10 19 is not attenuated 
substantially if sources are located at z < 0.03, and the 
flux at E = 5 x 10 19 is not attenuated if sources are 
located at z < 0.1. Assuming this range of redshifts 
for charged entries in Table I, we see that there is no 
problem to explain all of them except lines 3 and 5. The 
ray (21-22) has sufficiently low energy so that BL Lac 
B2 0804+35 with z=0.082 can be a source without any 
problem with attenuation. Likewise, the BL Lac 19 can 
be an actual source of the ray (19-20) if it has redshift 
z < 0.1. Thus, the correlations due to charged particles 
listed in Table I can be explained within the framework 
of conventional physics if sources with unknown redshifts 
are assumed to be within z < 0.1. 

Consider now the entries corresponding to neutral par- 
ticles. According to Ref. 0], photons can be UHECR 
primaries at any energy E > 10 19 eV provided that 
EGMF is smaller than 10 -9 G, energy spectrum at the 
source is hard, oc E~ a with a < 2 and the maximum en- 
ergy is high enough. (Note that correlations with charged 
particles imply that EGMF is in this range anyway.) 
Therefore, entries 3 and 17 can be explained by pho- 



tons. The real difficulty is with rays 5,7,8 and 16. Note, 
however, that according to MC simulations, there should 
be 7 background events in the Table I on average. 

To summarize, the idea of using the Galactic magnetic 
field as a mass-spectrograph of UHECR seems to work. 
The correlations between UHECR and BL Lacs substan- 
tially improve when arrival directions of cosmic rays are 
corrected for GMF. If not a statistical fluctuation, this 
implies the following: i) cosmic rays of highest energies 
contain a substantial fraction of protons ii) extragalac- 
tic magnetic fields have little effect on propagation of 
UHECR even from cosmological distances iii) the model 
of the Galactic magnetic field described by Eq. (||) is 
roughly correct. 

Finally, the significance of correlations with charged 
particles is p < 10 -3 . This may be considered as yet 
more evidence that BL Lacs are sources of UHECR. In- 
terestingly, the events and BL Lacs contributing into the 
correlation with lowest p found in Ref. || and the charged 
events of Table I do not overlap, so the two correlations 
should be considered as independent. 
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